!  Initial condition (density, magnetic field, velocity) 
!  for a particular configuration of magnetic carpets.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: linitial_condition = .true.
!
!***************************************************************
module InitialCondition
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet  
  use Mpicomm
  use Messages
  use Boundcond ! for the core boundary communication
!
  implicit none
!
  include '../initial_condition.h'
!
! ampl = amplitude of the magnetic field
!
  real :: ampl = 1.0
  character (len=labellen) :: config='parasitic_polarities'
!
  namelist /initial_condition_pars/ &
    ampl, config
!
  contains
!***********************************************************************
  subroutine register_initial_condition()
!
!  Configure pre-initialised (i.e. before parameter read) variables
!  which should be know to be able to evaluate
!
!  30-april-15/simon: coded
!
  endsubroutine register_initial_condition
!***********************************************************************
  subroutine initial_condition_uu(f)
!
!  Initialize the velocity field.
!
!  07-may-09/wlad: coded
!
    real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
    call keep_compiler_quiet(f)
!
  endsubroutine initial_condition_uu
!***********************************************************************
  subroutine initial_condition_lnrho(f)
!
!  Initialize logarithmic density. init_lnrho 
!  will take care of converting it to linear 
!  density if you use ldensity_nolog
!
!  07-may-09/wlad: coded
!
    real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
    call keep_compiler_quiet(f)
!
  endsubroutine initial_condition_lnrho
!***********************************************************************
  subroutine initial_condition_aa(f)
!
!  Initialize the magnetic vector potential.
!
!  30-april-15/simon: coded
!
!  Magnetic carpet consisting of three spots.
!
!  Created 2015-04-30 by Simon Candelaresi (Iomsn)
!
    use Mpicomm, only: stop_it
    use Poisson
    use Sub
!    
    real, dimension (mx,my,mz,mfarray) :: f
!        
!   The next variables are used for the uncurling.
    integer :: l, j, ju, k
    real, dimension (nx,ny,nz,3) :: jj, tmpJ  ! This is phi for poisson.f90
!
!   clear the magnetic field to zero
    f(:,:,:,iax:iaz) = 0.
!
    if (config == 'homogeneous') then
        do m = 1, my, 1
            do l = 1, mx, 1
                f(l,m,:,iax) = -ampl*y(m)/2
                f(l,m,:,iay) = ampl*x(l)/2
            enddo
        enddo
    endif
    
    if (config == 'parasitic_polarities') then
        do n = 1, mz, 1
            do m = 1, my, 1
                do l = 1, mx, 1
                    f(l,m,n,iax) = ampl * (&
                        2. * ((x(l) - 10.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2)&
                        ** (-3.0 / 2.) * y(m) + 2. * ((x(l) + 8.) ** 2 + y(m) **&
                        2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * y(m) + 2. * ((x(l)&
                        + 10.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
                        / 2.) * y(m) + 2. * ((x(l) + 6.) ** 2 + y(m) ** 2 + (z(n)&
                        + 0.85) ** 2) ** (-3.0 / 2.) * y(m) + 2. * (x(l) ** 2 + y(m)&
                        ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * y(m) + 2. *&
                        ((x(l) + 2.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
                        / 2.) * y(m) + 2. * ((x(l) - 2.) ** 2 + y(m) ** 2 + (z(n)&
                        + 0.85) ** 2) ** (-3.0 / 2.) * y(m) + 2. * ((x(l) - 8.)&
                        ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 /&
                        2.) * (y(m) - 8.) + 2. * ((x(l) - 6.) ** 2 + (y(m) - 8.)&
                        ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) - 8.) +&
                        2. * ((x(l) - 10.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85)&
                        ** 2) ** (-3.0 / 2.) * (y(m) - 8.) + 2. * ((x(l) - 8.) **&
                        2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.)&
                        * (y(m) + 8.) + 2. * ((x(l) - 6.) ** 2 + (y(m) + 8.) ** 2&
                        + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) + 8.) + 2. *&
                        ((x(l) - 10.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) **&
                        2) ** (-3.0 / 2.) * (y(m) + 8.) + 2. * ((x(l) + 8.) ** 2 +&
                        (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m)&
                        - 8.) + 2. * ((x(l) + 10.) ** 2 + (y(m) - 8.) ** 2 +&
                        (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) - 8.) + 2. * ((x(l)&
                        + 6.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) **&
                        (-3.0 / 2.) * (y(m) - 8.) + 2. * ((x(l) + 8.) ** 2 + (y(m)&
                        + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m)&
                        + 8.) + 2. * ((x(l) + 10.) ** 2 + (y(m) + 8.) ** 2 + (z(n)&
                        + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) + 8.) + 2. * ((x(l)&
                        + 6.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
                        / 2.) * (y(m) + 8.) + 2. * (x(l) ** 2 + (y(m) + 8.) **&
                        2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) + 8.) + 2.&
                        * ((x(l) + 2.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) **&
                        2) ** (-3.0 / 2.) * (y(m) + 8.) + 2. * ((x(l) - 2.) ** 2&
                        + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) *&
                        (y(m) + 8.) + 2. * (x(l) ** 2 + (y(m) - 8.) ** 2 + (z(n) +&
                        0.85) ** 2) ** (-3.0 / 2.) * (y(m) - 8.) + 2. * ((x(l) +&
                        2.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
                        / 2.) * (y(m) - 8.) + 2. * ((x(l) - 2.) ** 2 + (y(m) - 8.)&
                        ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) - 8.)&
                        + 2. * ((x(l) - 8.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2)&
                        ** (-3.0 / 2.) * y(m) + 2. * ((x(l) - 6.) ** 2 + y(m) **&
                        2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * y(m))
                    f(l,m,n,iay) = ampl * (&
                        -2. * ((x(l) - 10.) ** 2 + y(m) ** 2 + (z(n) + 0.85) **&
                        2) ** (-3.0 / 2.) * (x(l) - 10.) - 2. * ((x(l) + 8.) ** 2&
                        + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) +&
                        8.) - 2. * ((x(l) + 10.) ** 2 + y(m) ** 2 + (z(n) + 0.85)&
                        ** 2) ** (-3.0 / 2.) * (x(l) + 10.) - 2. * ((x(l) + 6.) **&
                        2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l)&
                        + 6.) - 2. * (x(l) ** 2 + y(m) ** 2 + (z(n) + 0.85) **&
                        2) ** (-3.0 / 2.) * x(l) - 2. * ((x(l) + 2.) ** 2 + y(m) **&
                        2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 2.) - 2.&
                        * ((x(l) - 2.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) **&
                        (-3.0 / 2.) * (x(l) - 2.) - 2. * ((x(l) - 8.) ** 2 + (y(m)&
                        - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l)&
                        - 8.) - 2. * ((x(l) - 6.) ** 2 + (y(m) - 8.) ** 2 + (z(n)&
                        + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) - 6.) - 2. * ((x(l) -&
                        10.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
                        / 2.) * (x(l) - 10.) - 2. * ((x(l) - 8.) ** 2 + (y(m) +&
                        8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) - 8.)&
                        - 2. * ((x(l) - 6.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85)&
                        ** 2) ** (-3.0 / 2.) * (x(l) - 6.) - 2. * ((x(l) - 10.)&
                        ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
                        / 2.) * (x(l) - 10.) - 2. * ((x(l) + 8.) ** 2 + (y(m) - 8.)&
                        ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 8.)&
                        - 2. * ((x(l) + 10.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85)&
                        ** 2) ** (-3.0 / 2.) * (x(l) + 10.) - 2. * ((x(l) + 6.)&
                        ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 /&
                        2.) * (x(l) + 6.) - 2. * ((x(l) + 8.) ** 2 + (y(m) + 8.) **&
                        2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) - 2.&
                        * ((x(l) + 10.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85)&
                        ** 2) ** (-3.0 / 2.) * (x(l) + 10.) - 2. * ((x(l) + 6.) **&
                        2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.)&
                        * (x(l) + 6.) - 2. * (x(l) ** 2 + (y(m) + 8.) ** 2 + (z(n)&
                        + 0.85) ** 2) ** (-3.0 / 2.) * x(l) - 2. * ((x(l) + 2.)&
                        ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.)&
                        * (x(l) + 2.) - 2. * ((x(l) - 2.) ** 2 + (y(m) + 8.) **&
                        2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) - 2.) - 2.&
                        * (x(l) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) **&
                        (-3.0 / 2.) * x(l) - 2. * ((x(l) + 2.) ** 2 + (y(m) - 8.)&
                        ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 2.) -&
                        2. * ((x(l) - 2.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85)&
                        ** 2) ** (-3.0 / 2.) * (x(l) - 2.) - 2. * ((x(l) - 8.) **&
                        2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l)&
                        - 8.) - 2. * ((x(l) - 6.) ** 2 + y(m) ** 2 + (z(n) + 0.85)&
                        ** 2) ** (-3.0 / 2.) * (x(l) - 6.) + x(l))
                    f(l,m,n,iaz) = 0
                enddo          
            enddo          
        enddo

    elseif (config == 'single_parasitic_polarity') then
        do n = 1, mz, 1
            do m = 1, my, 1
                do l = 1, mx, 1
                    f(l,m,n,iax) = ampl * (&        
                        2. * (x(l) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
                        / 2.) * y(m) + 2. * ((x(l) + 8.) ** 2 + y(m) ** 2 + (z(n)&
                        + 0.85) ** 2) ** (-3.0 / 2.) * y(m) + 2. * ((x(l) - 8.)&
                        ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * y(m)&
                        + 2. * (x(l) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) **&
                        2) ** (-3.0 / 2.) * (y(m) - 8.) + 2. * (x(l) ** 2 + (y(m)&
                        + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m)&
                        + 8.) + 2. * ((x(l) + 8.) ** 2 + (y(m) + 8.) ** 2 + (z(n)&
                        + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) + 8.) + 2. * ((x(l) +&
                        8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
                        / 2.) * (y(m) - 8.) + 2. * ((x(l) - 8.) ** 2 + (y(m) + 8.)&
                        ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) + 8.)&
                        + 2. * ((x(l) - 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85)&
                        ** 2) ** (-3.0 / 2.) * (y(m) - 8.))
                    f(l,m,n,iay) = ampl * (&
                        x(l) - 2. * (x(l) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2)&
                        ** (-3.0 / 2.) * x(l) - 2. * ((x(l) + 8.) ** 2 + y(m) **&
                        2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) - 2.&
                        * ((x(l) - 8.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
                        / 2.) * (x(l) - 8.) - 2. * (x(l) ** 2 + (y(m) - 8.) **&
                        2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * x(l) - 2. * (x(l)&
                        ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
                        / 2.) * x(l) - 2. * ((x(l) + 8.) ** 2 + (y(m) + 8.) ** 2 +&
                        (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) - 2. * ((x(l)&
                        + 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2)&
                        ** (-3.0 / 2.) * (x(l) + 8.) - 2. * ((x(l) - 8.) ** 2 + (y(m)&
                        + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l)&
                        - 8.) - 2. * ((x(l) - 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n)&
                        + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) - 8.))
                    f(l,m,n,iaz) = 0
                enddo          
            enddo          
        enddo

    elseif (config == 'dominant_polarities') then
        do n = 1, mz, 1
            do m = 1, my, 1
                do l = 1, mx, 1
                    f(l,m,n,iax) = ampl * (&
                        -0.3 * ((x(l) - 105) ** 2 + ((y(m) + 8) ** 2) + (z(n) +&
                        0.3) ** 2) ** (-3.0 / 2.) * (y(m) + 8) - 0.3 * ((x(l) - 8.)&
                        ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 /&
                        2.) * (y(m) + 8) - 0.3 * ((x(l) + 8.) ** 2 + ((y(m) - 8)&
                        ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) - 8) - 0.3&
                        * ((x(l) + 105) ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3)&
                        ** 2) ** (-3.0 / 2.) * (y(m) - 8) - 0.3 * ((x(l) + 5.5) **&
                        2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.)&
                        * (y(m) - 8) - 0.3 * ((x(l) + 8.) ** 2 + ((y(m) + 8) **&
                        2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) + 8) - 0.3&
                        * ((x(l) + 105) ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3) **&
                        2) ** (-3.0 / 2.) * (y(m) + 8) - 0.3 * ((x(l) + 5.5) ** 2&
                        + ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) *&
                        (y(m) + 8) - 0.3 * (x(l) ** 2 + ((y(m) + 8) ** 2) + (z(n)&
                        + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) + 8) - 0.3 * ((x(l) +&
                        2.5) ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
                        / 2.) * (y(m) + 8) - 0.3 * ((x(l) - 2.5) ** 2 + ((y(m)&
                        + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) + 8)&
                        - 0.3 * (x(l) ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) **&
                        2) ** (-3.0 / 2.) * (y(m) - 8) - 0.3 * ((x(l) + 2.5) ** 2&
                        + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) *&
                        (y(m) - 8) - 0.3 * ((x(l) - 2.5) ** 2 + ((y(m) - 8) ** 2)&
                        + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) - 8) - 0.3 *&
                        ((x(l) - 8.) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
                        / 2.) * (y(m)) - 0.3 * ((x(l) - 5.5) ** 2 + (y(m) ** 2)&
                        + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m)) - 0.3 * ((x(l)&
                        - 105) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
                        / 2.) * (y(m)) - 0.3 * ((x(l) + 8.) ** 2 + (y(m) ** 2) +&
                        (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m)) - 0.3 * ((x(l)&
                        + 105) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 /&
                        2.) * (y(m)) - 0.3 * ((x(l) + 5.5) ** 2 + (y(m) ** 2) + (z(n)&
                        + 0.3) ** 2) ** (-3.0 / 2.) * (y(m)) - 0.3 * ((x(l) -&
                        2.5) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.)&
                        * (y(m)) - 0.3 * (x(l) ** 2 + (y(m) ** 2) + (z(n) + 0.3)&
                        ** 2) ** (-3.0 / 2.) * (y(m)) - 0.3 * ((x(l) + 2.5) ** 2&
                        + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m))&
                        - 0.3 * ((x(l) - 8.) ** 2 + ((y(m) - 8) ** 2) + (z(n) +&
                        0.3) ** 2) ** (-3.0 / 2.) * (y(m) - 8) - 0.3 * ((x(l) - 5.5)&
                        ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
                        / 2.) * (y(m) - 8) - 0.3 * ((x(l) - 105) ** 2 + ((y(m) - 8)&
                        ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) - 8) -&
                        0.3 * ((x(l) - 5.5) ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3)&
                        ** 2) ** (-3.0 / 2.) * (y(m) + 8))
                    f(l,m,n,iay) = ampl * (&
                        0.3 * ((x(l) - 8.) ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3)&
                        ** 2) ** (-3.0 / 2.) * (x(l) - 8.) + 0.3 * ((x(l) - 5.5)&
                        ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 /&
                        2.) * (x(l) - 5.5) + 0.3 * ((x(l) - 105) ** 2 + ((y(m) -&
                        8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 105)&
                        + 0.3 * ((x(l) - 8.) ** 2 + ((y(m) + 8) ** 2) + (z(n) +&
                        0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 8.) + 0.3 * ((x(l) - 5.5)&
                        ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
                        / 2.) * (x(l) - 5.5) + 0.3 * ((x(l) - 105) ** 2 + ((y(m)&
                        + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 105)&
                        + 0.3 * ((x(l) + 8.) ** 2 + ((y(m) - 8) ** 2) + (z(n)&
                        + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) + 0.3 * ((x(l) +&
                        105) ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
                        / 2.) * (x(l) + 105) + 0.3 * ((x(l) + 5.5) ** 2 + ((y(m)&
                        - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) +&
                        5.5) + 0.3 * ((x(l) + 8.) ** 2 + ((y(m) + 8) ** 2) + (z(n)&
                        + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) + 0.3 * ((x(l)&
                        + 105) ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
                        / 2.) * (x(l) + 105) + 0.3 * ((x(l) + 5.5) ** 2 + ((y(m)&
                        + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l)&
                        + 5.5) + 0.3 * (x(l) ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3)&
                        ** 2) ** (-3.0 / 2.) * x(l) + 0.3 * ((x(l) + 2.5) ** 2&
                        + ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) *&
                        (x(l) + 2.5) + 0.3 * ((x(l) - 2.5) ** 2 + ((y(m) + 8) **&
                        2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 2.5) + 0.5D-1&
                        * x(l) + 0.3 * (x(l) ** 2 + ((y(m) - 8) ** 2) + (z(n)&
                        + 0.3) ** 2) ** (-3.0 / 2.) * x(l) + 0.3 * ((x(l) + 2.5)&
                        ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.)&
                        * (x(l) + 2.5) + 0.3 * ((x(l) - 2.5) ** 2 + ((y(m) - 8)&
                        ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 2.5)&
                        + 0.3 * ((x(l) - 5.5) ** 2 + (y(m) ** 2) + (z(n) + 0.3) **&
                        2) ** (-3.0 / 2.) * (x(l) - 5.5) + 0.3 * ((x(l) - 105) **&
                        2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l)&
                        - 105) + 0.3 * ((x(l) - 8.) ** 2 + (y(m) ** 2) + (z(n)&
                        + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 8.) + 0.3 * ((x(l) +&
                        8.) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.)&
                        * (x(l) + 8.) + 0.3 * ((x(l) + 105) ** 2 + (y(m) ** 2) +&
                        (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) + 105) + 0.3 *&
                        ((x(l) + 5.5) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
                        / 2.) * (x(l) + 5.5) + 0.3 * (x(l) ** 2 + (y(m) ** 2)&
                        + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * x(l) + 0.3 * ((x(l)&
                        + 2.5) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 /&
                        2.) * (x(l) + 2.5) + 0.3 * ((x(l) - 2.5) ** 2 + (y(m) **&
                        2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 2.5))
                    f(l,m,n,iaz) = 0
                enddo          
            enddo          
        enddo

    elseif (config == 'single_dominant_polarity') then
        do n = 1, mz, 1
            do m = 1, my, 1
                do l = 1, mx, 1
                    f(l,m,n,iax) = ampl * (&        
                        -1. * (x(l) ** 2 + y(m) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0&
                        / 2.) * y(m) - 1. * ((x(l) + 8.) ** 2 + y(m) ** 2 + (z(n)&
                        + 0.5) ** 2) ** (-3.0 / 2.) * y(m) - 1. * ((x(l) - 8.)&
                        ** 2 + y(m) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * y(m)&
                        - 1. * (x(l) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.5) ** 2)&
                        ** (-3.0 / 2.) * (y(m) - 8.) - 1. * (x(l) ** 2 + (y(m) +&
                        8.) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * (y(m) + 8.)&
                        - 1. * ((x(l) + 8.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.5)&
                        ** 2) ** (-3.0 / 2.) * (y(m) + 8.) - 1. * ((x(l) + 8.)&
                        ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.)&
                        * (y(m) - 8.) - 1. * ((x(l) - 8.) ** 2 + (y(m) + 8.) **&
                        2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * (y(m) + 8.) - 1. *&
                        ((x(l) - 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.5) ** 2)&
                        ** (-3.0 / 2.) * (y(m) - 8.))
                    f(l,m,n,iay) = ampl * (&
                        x(l) + 1. * (x(l) ** 2 + y(m) ** 2 + (z(n) + 0.5) ** 2)&
                        ** (-3.0 / 2.) * x(l) + 1. * ((x(l) + 8.) ** 2 + y(m) ** 2&
                        + (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) + 1. *&
                        ((x(l) - 8.) ** 2 + y(m) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0&
                        / 2.) * (x(l) - 8.) + 1. * (x(l) ** 2 + (y(m) - 8.) ** 2&
                        + (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * x(l) + 1. * (x(l) **&
                        2 + (y(m) + 8.) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.)&
                        * x(l) + 1. * ((x(l) + 8.) ** 2 + (y(m) + 8.) ** 2 + (z(n)&
                        + 0.5) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) + 1. * ((x(l)&
                        + 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0&
                        / 2.) * (x(l) + 8.) + 1. * ((x(l) - 8.) ** 2 + (y(m) + 8.)&
                        ** 2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * (x(l) - 8.)&
                        + 1. * ((x(l) - 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.5)&
                        ** 2) ** (-3.0 / 2.) * (x(l) - 8.))
                    f(l,m,n,iaz) = 0
                enddo          
            enddo          
        enddo

    endif
!
  endsubroutine initial_condition_aa
!***********************************************************************
    subroutine read_initial_condition_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=initial_condition_pars, IOSTAT=iostat)
!
    endsubroutine read_initial_condition_pars
!***********************************************************************
    subroutine write_initial_condition_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=initial_condition_pars)
!
    endsubroutine write_initial_condition_pars
!***********************************************************************
!***********************************************************************
!************        DO NOT DELETE THE FOLLOWING       **************
!********************************************************************
!**  This is an automatically generated include file that creates  **
!**  copies dummy routines from noinitial_condition.f90 for any    **
!**  InitialCondition routines not implemented in this file        **
!**                                                                **
  include '../initial_condition_dummies.inc'
!********************************************************************
endmodule InitialCondition
